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Abstract 

One of the few methods for generating efficient function spaces for multi-D Schrodinger eigen- 
problems is given by Garashchuk and Light in J.Chem.Phys. 114 (2001) 3929. Their Gaussian 
basis functions are wider and sparser in high potential regions, and narrower and denser in low 
ones. We suggest a modification of their approach based on the following observation: In very 
steep potential regions, wide, sparse, Gaussians should be avoided even if their centers have high 
potential values. Our numerical results illustrate that a dramatic improvement in accuracy may be 
obtained in this way. We also compare the errors of collocation to those of a Galerkin approach, 
test a criterion for scaling Gaussian widths based on deviation from orthogonality of collocation 
eigenfunctions, and suggest a criterion for scaling Gaussian widths based on Hamiltonian trace 
minimization. 



1 Introduction 

Numerical prediction of the properties and behaviour of quantum mechanical systems is a long standing 
challenge. A major difficulty stems from the high dimensions of such systems. An n > 3 particle 
mechanical system generally has 3n — 6 degrees of freedom, and so to calculate the state of an n > 3 
particle quantum system we must solve a PDE in 3n— 6 variables. Presently the limit is n = 5 particles 
{i.e. 9 variables) in the largest systems for which many eigenvalues/eigenstates were obtained from 
Schrodinger's equation pp. However, there are many larger systems for which numerical predictions 
would be very useful for chemistry biology (where many important molecules have more than five 
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atoms), and for the emerging field of nano engineering. Clearly, we would like to push the limit 
further. 

To efficiently treat multi-D (i.e. multi variable) PDEs we must consider the function spaces where 
approximations are constructed. Naively, we can form such function spaces by taking tensor products 
of 1-D function spaces. Many of the 1-D spaces standardly used (e.g. Fourier spaces, spaces of 
orthogonal polynomials of different kinds) have a basis of localized functions peaked above the nodes 
of a 1-D quadrature formula, these are called 1-D DVRs [2], [3]. The corresponding multi-D DVR 
basis of the product space consists of localized multi variable functions peaked above the points of Q, 
the (hyper) rectangular array of nodes of the associated product cubature formula. Typically, at least 
10 1-D DVR basis functions are taken for each variable, so the dimension of the product space is at 
least 10 3n ~ 6 . This gets much too large, even when the number of particles n is not too big. 

DVR bases can give us some insight regarding the inefficiency of multi-D tensor product spaces. 
Bound solutions of Schrodinger's eigenproblem decay exponentially in high potential regions. Thus, 
for a given system we can apriori estimate f2, the region where wavefunctions up to a given energy 
(eigenvalue) are supported. Typically, f2 is non rectangular and covering it with the rectangular grid 
Q is inefficient. Many nodes will be outside f2 and therefore many of the multi-D product DVR basis 
functions will contribute little to the approximation of our wavefunctions. Evidently, efficient localized 
basis functions should be peaked in Q, and not outside it. Several approaches have been suggested for 
obtaining such basis functions. 

One very useful solution is to prune a product DVR basis by discarding functions whose nodes lie 
in high potential regions [3]. Another approach, still in its initial stages, is to construct DVRs based on 
non product cubatures with various non rectangular node distributions [I]. In both cases the resulting 
DVR basis functions are orthogonal. However, in these approaches it is not easy to refine the node 
distribution within Q according to predictable wavefunction properties. Generally, wavefunctions 
are more oscillatory in low potential regions compared to high ones and node densities in $7 should 
vary accordingly. Therefore we may try to construct function spaces by appropriately distributing 
nodes and then placing a localized basis function above each. Gaussians [3] [U] [7J, B-splines [8], and 
inverse multi quadrics [11], have been used in this way. Although the resulting basis functions are not 
orthogonal, the flexible choice of nodes and basis function shape parameters can give highly accurate 
results. This has been illustrated in 1-D problems [5] [6]. However the only general method we are 
aware of for constructing multi-D bases of localized functions whose node density varies with potential 
values is given by Garashchuk and Light [7J. 

Here we closely examine the method of [7j for choosing nodes and widths of Gaussian basis func- 
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tions. Their basis functions are narrow and close together in low potential regions, while they are 
wide and sparsely distributed in high potential regions. We believe that this needs to be refined in 
steep potential regions, where points with low and high potential values are very close together. This 
situation is common in molecular potential surfaces for small nuclear separation. We suggest a simple 
modification of the nodes/widths choice of [7J. Its essence is to avoid overly wide and sparse basis 
functions in hard potential regions. We illustrate that this can give a dramatic improvement in accu- 
racy by calculating the eigenfucntions and eigenvalues of two 1-D Morse Hamiltonians. Our results 
indicate that the algorithm of [7] should be adapted along these lines in the multi-D case. 

Our work also explores the use of collocation with Gaussian basis functions. Using collocation, 
rather than the standard Galerkin approach, is very appealing since calculating integrals of potential 
matrix elements is avoided. Yang and Peet have previously used collocation with Gaussian basis 
functions [13], and in their calculations performance was similiar to the Galerkin approach. They have 
used the same basis functions with both methods. However in our calculations basis function widths 
were adjusted optimally for each method and Galerkin accuracy was much better than collocation. 
Collocation has many advantages and it can give high quality results as shown in the work of Yang 
and Peet [13], [H], [15], and in many papers based on it. However, our results illustrate that caution 
should be exercised before embracing collocation as equally accurate to Galerkin. 

Other aspects we examine are the use of collocation with the quasi randomly distributed Gaussians 
of [TJ, and using the deviation from orthogonality of eigenfunctions calculated by collocation to guide 
us in the choice of basis function widths. 

The content of this paper follows. In section [2] we describe the distributed Gaussian bases used 
for Schrodinger's equation leading up to [7J, then we overview the collocation approach. In section 
[3] we discuss the nodes/widths choice of [7J, suggest improvements, and raise a few further questions 
about collocation and Gaussian basis functions. Section H] gives numerical results for two 1-D Morse 
Hamiltonians. Finally, we conclude and suggest some future directions in section [5j 

2 Distributed Gaussian bases and collocation 

Here we discuss the application of distributed Gaussian bases and collocation in the Schrodinger 
eigenproblem, i.e. the problem of finding eigenvalues and eigenfunctions of the Hamiltonian operator 
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We concentrate on the problem of finding all bound states below an energy cutoff E max , i.e. all 
eigenfunctions {V'n} and eigenvalues {E n } belonging to the discrete part of the spectrum with E n < 

E m ax ■ 

Several researchers have studied the use of Gaussian basis functions for this proble We pick up 
an important thread in this research by considering the work of Hamilton and Light [5 J , and continuing 
to Poirier and Light [B|, and to Garashchuk and Light [7j. 

The basis functions used in [5] are real Gaussians 

< Pi (x)=r H e- cw *<- x - x 4' i = l,...,n (2) 
whose parameters are named as follows. 

• Xi is called the center or node of ipi. 

• u>i > is the width parameter of ipi. 

• c > is the "global" width parameter. 

• The normalization parameter rji is chosen so that {tpi\<Pi) = 1. 

Much depends, of course, on the choice of the In + 1 parameters c, {wi}, {x^}. For 1-D problems 
classical mechanics arguments are used in [5] to construct very efficient bases, "tailored" to a given 
problem. However, such arguments are not readily extended to multi-D, and a simpler approach is 
applied in [5] in this case. Its rationale can be explained as follows. The highest eigenfunction with 
energy (approximately) E max is the most oscillatory in low potential regions, while reaching furthest 
out to high potential regions. Therefore a basis suitable for its approximation will be suitable for 
all lower eigenfunctions. The most severe oscillations appear near the potential minimum where the 
local De-Brogue wavelength is 2ir / \J 2m{E max — V m i n )/fi 2 . Thus the {x{} are chosen in |i5j as points 

1 A thorough discussion of other types of basis functions is beyond the scope of this work. However, the performance of 
well chosen Gaussian bases is good. In [8] Galerkin b-spline calculations for Morse Hamiltonians are reported. Generally, 
results obtained there with approximately 100 B-splines have much greater errors compared to results obtained with 
48 Gaussians in [5], [7], and in section [4] here (to be fair we note that the Morse Hamiltonians were not identical). 
Moreover, it is illustrated numerically in the first reference of [8] that eigenvalue errors using n Gaussians scale as a~ n 
while the error using n B-splines of order k scales as (^) 2fc , i.e. the Gaussian error converges asymptotically faster. 
A thorough analytical and numerical investigation of using high order B-splines in Schrodingers equation is certainly 
in place. However, we presently opt for Gaussians because of the reasons given above, and because the lessons on 
node/ width distributions we obtain using Gaussians will probably be useful for other types of basis functions including 
B-splines. 
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of a regular grid whose spacing is at most 1/4-7T of this De-Broglie wavelength (any sufficiently fine 
spacing would do). On the other hand, the highest eigenfunction is exponentially decaying in the 
classically forbidden region where V(x) > E max , thus the X{ are confined to the classically accessible 
region V(x) < E max . Regular grids are used in [5] which allow an easy notion of grid spacing, however 
this is by no means a restriction to rectangular grids. For the Henon-Heiles problem the Xi chosen 
in [5] are points of a regular hexagonal grid pruned to fit within the classically accessible triangular 
region. Regarding width parameters, the choice W{ = 1/Ax 2 , is advocated in [5]. It is interesting to 
note that this choice arises if we require the curvature of the cpi at their centers x\ (i.e. ^^\ x=x .) to 
be compatible with the maximal curvature of the highest eigenfunction. Having chosen the Wi we can 
still vary the global width parameter c. Let us define the matrix S by 



In [5] the most accurate Hamiltonian eigenvalues are obtained for low c values for which the condition 
number of S is large while still allowing stable numerical calculations. 

Hamilton and Light [5] provide efficient Gaussian bases whose centers occupy only relevant regions 
of configuration space. However, for multi-D problems the density of centers in [5] is constant, and 
dictated by the highest oscillations in the problem. Such high center density is necessary only in the 
lowest potential regions. As in the 1-D examples in [5] we expect that sparser Gaussian centers can 
be taken in higher potential regions. This is the subject of [TJ, but we first discuss [6J which provides 
the background for [TJ, and which unveils important, and surprising, properties of Gaussian bases. 

Intuitively, we tend to reject badly conditioned basis functions which are "almost linearly depen- 
dent". Thus we may expect that "good" Gaussian basis functions will not be "too wide" and their 
centers will not be "too close" together. In [8J Poirier and Light show the contrary; that "good" 
Gaussian bases can actually have clustered centers and high condition numbers. They regard the 
choice of Gaussian basis parameters as an optimization problem. Let us define the matrix H by 



Poirier and Light solve for widths and centers which minimize the functional trace(S~ l H). This is 
based on the "variational principle" : each eigenvalue of H is smaller than the corresponding eigenvalue 
of see MacDonald [9j. Therefore minimizing trace(S~ l H) is equivalent to minimizing eigenvalue 

errors. Note that S^ 1 !! represents H as a linear operator in the non-orthogonal Gaussian basis, while 
H represents the bilinear form (-\H-). Note also that including S^ 1 in our functional we ensure that the 
Gaussians yielded by the optimization will be linearly independent, however their condition number 
is not constrained. 




(3) 




(4) 
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An illumating example is provided by the 1-D harmonic oscillator, H = —^-§^1 + ^muj 2 x 2 . 
Rather than giving a "well conditioned" basis, the optimization produces Gaussians which are small 

^ 2 

perturbations of the harmonic oscillator ground state ipo(x) = k$e~~ x , where a = \frnu)Jh. These 
basis functions have large overlaps, and it is reported in [6] that they reproduce the harmonic oscillator 

eigenvalues to machine accuracy. At first this may seem surprising, but a simple explanation is 

1-1 _«i 2 

provided in [6j. The n'th harmonic oscillator eigenstate is ip n (x) = k n e 2 x H n (ax), where H n is the 

nth Hermite polynomial. Using the formula e~ y2 H n (y) = (— l)™^7re _s/2 (from e.g. flO]) we see that 

jn oP 2 / \ ^ cP 2 jfc oP 2 

Jj^—Q-— x = e~~ x H n (^x), and therefore ip n is in span{-^jre~~ x \k = 0, ... , n}. Consider 

now the space mentioned above of Gaussians closely clustered near the potential minimum. It contains 

flk cP 2 

difference approximations to the derivatives - d ~^^~~ x > an d since the centers are close together the 
accuracy of these difference approximations is high. Therefore n closely clustered translations of the 
ground state span a linear space which gives excellent approximations of the first n harmonic oscillator 
eigenstates. Clustering of optimal Gaussian centers is observed in [B] also for a Morse potential, 
however in this case there are several clustering locations away from the potential minimum. The 
high condition number of such optimal Gaussian bases can be problematic in practical large scale 
calculations. However, the harmonic oscillator example indicates that high overlaps are needed only 
to produce accurate approximations of Gaussian derivatives. Poirier and Light therefore propose to 
include derivatives of Gaussians among our basis functions to begin with, and so avoid the need for 
large Gaussian overlaps. This direction is not pursued further in [6]. Instead, the bad conditioning of 
optimal Gaussian bases is addressed by introducing a modified functional, 

tracers- 1 H) + A||5- J|| . (5) 

Gaussians optimizing this functional strike a compromise, depending on the value of A, between 
accuracy and conditioning. A similiar functional provides the starting point for [7]. 

In [7] Garashchuk and Light suggest a practical way of constructing multi-D Gaussian bases, whose 
node density varies with potential values. They define the functional 

trace(H)+\J2T%- , (6) 

which is different from the functional © used in [6]. Note particularly that trace{S~ 1 H) from ([5]) 
is replaced by trace(H) in ([6]). This is minimized for 1-D Morse and double well Hamiltonians to 
find the density of optimal Gaussian centers p(xi) = 2(x- + i-x- 1) anc ^ th e fr widths Wi. It is then 
observed that p{xi) and Wi fit linear functions of the potential, i.e. there are constants ao,ax,bo,bi 
such that p{xi) ~ ao(ai — V{xi)), and Wi ~ &o(&l — U(xj)). A similiar optimization for realistic large 
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scale problems would be very time consuming. However, Garashchuk and Light regard these two 1-D 
examples as manifestations of a general rule: "Gaussian bases which have center density and width 
parameters fitting appropriate linear functions of the potential will be good". It is also observed 
in [7J that Hamiltonian eigenvalue errors decrease when the global width parameter c is decreased 
until unstable numerical calculations set in when cond(S) > 10 12 . Based on these observations a 
simple procedure for choosing multi-D Gaussian basis functions is suggested in [7J. Fixing E max and 
a parameter p they define the density function 

Kx) - f *-- y W + " V . (7) 



With 7 = 1 the density p depends linearly on the potential. The value 7 = 1/2 is also considered 
in [TJ, in this case p corresponds to the local De-Broglie wave length of the highest eigenf unction. 
Having defined p(x) candidate points x are quasi-randomly generated in the region V(x) < E max . x 
is accepted as a Gaussian center if 

p(x) > r (8) 

where r G [0, 1] is random. With the centers at hand, the widths are defined 

Wi = E max - V(xi) + p > . (9) 

Then c is taken as small as possible until cond(S) ~ 10 12 . 

This approach is applied in [7J to eigenvalue calculations of several multi-D Hamiltonians. The 
results are much better, or not worse, compared to those of previous methods. Therefore the quasi 
random distributed Gaussians of [7J appear to be a very good choice. It is our purpose here to make 
some further observations and suggest a few possible refinements of the approach of [7J. Among other 
things we shall use the Gaussian basis functions of [7J in a collocation algorithm, which we discuss 
next. 

The collocation approach for finding the eigenvalues and eigenfunctions of a Hamiltonian H may be 
described as follows. First, we choose basis functions tpi, . . . , <p n , and a set of points x%, . . . ,x n called 
collocation nodes. We then seek functions ip = J^iLi^iVi) an d eigenvalues E, such that (Hip)(xi) = 
Eifj(xi), i = 1, . . . ,n. This gives the following generalized matrix eigenproblem 

( -— A$ + y$ ) u = E®u, (10) 
V 2m J 

where $y = (fj(xi), (A0)y = Aifj(xi), V{j = 5ijV(xi). Since the exact eigenvalues are real we reject 
solutions with complex E, if there are any. Assuming <3? is invertible (|10p is equivalent to 

2m 



h 

$ -1 (A$) + $- 1 V$)u = Eu (11) 
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and changing coordinates to y = &u we obtain 



( 



(12) 



The matrix on the l.h.s. of (jlip is the collocation representation of the Hamiltonian in the ipi basis, 
while the matrix on the l.h.s. of (]12|) is the collocation representation in a basis 0j, spanning the same 
function space, and consisting of functions which are 1 on "their" node and zero on all other nodes, 



The collocation approach for spatial discretization of differential operators has a long history |16j . 
An important mile stone was the introduction of radial basis function (RBF) collocation by Kansa 
|17] . an approach which found much use in the PDEs of classical physics. For recent reviews on radial 
basis functions and their applications see [19]. RBFs have the form cpi(x,Wi) = f(\\x — Xi\\,Wi), where 

2 

/ is a real valued function peaked at and depending on the parameter vector w, e.g. f(r, w) = e~ wr 
giving Gaussian RBFs, or f(r, [w a ,wp]) = (1 + w 2 t r 2 )~ w i 3 ^ 2 giving generalized inverse multi quadric 
(GIMQ) RBFs (here wp > d where x € JR rf ). The centers xi are usually chosen as the collocation 
nodes. The attraction of RBF collocation stems from its great flexibility and simplicity. The basis 
functions shape and location can be adapted to a given problem by carefully choosing parameters. 
Moreover, the collocation equations are very simple to set up, they do not require calculation of multi 
dimensional integrals. 

The need for efficient methods is particularly acute for high dimensional problems in quantum 
mechanics. Yang and Peet were the first to show that collocation is feasible for the Schrodinger 
eigenproblem |13j . They applied collocation to 1-D Morse Hamiltonians together with the Gaussian 
RBFs of [5], and found that results were as good as the Galerkin approach with the same basis func- 
tions. The RBF collocation of [T_2] has become a standard method for 1-D Schrodinger eigenproblems 
(appearing as components in more complicated calculations). In multi-D problems collocation did 
not go, to the best of our knowledge, beyond the stage of initial explorations. Collocation was used 
for a 2-D problem (Ar-HCl Hamiltonian) in p3] with rectangular grids of nodes and tensor product 
function spaces. Again, performance was similiar to that of a Galerkin approach with the same basis 
functions. In [T5] collocation was used for the same 2-D problem with Gaussian RBFs whose centers 
were obtained by truncation of a rectangular grid outside the classically accessible region. However, in 
|15] only the first three eigenvalues were calculated using an iterative scheme based on projections to 
low dimensional subspaces. The only example we are aware of for calculation of high (additionally to 
low) lying eigenvalues and eigenfunctions of multi-D Hamiltonians using RBF collocation is given by 
Hu, Ho, and Rabitz [IT]. The nodes used in [JTJ were on hexagonal and triangular grids and the basis 
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functions were the GIMQ RBFs mentioned above. No direct comparison with the Galerkin approach 
was given in 

Having surveyed the topics of Gaussian RBFs and collocation for the Schrodinger eigenproblem 
we can now venture to take a closer look at some aspects of these topics. 

3 Observations and suggestions 

Here we give some observations and suggestions which arose while studying the application of Gaussian 
RBF collocation to Schrodinger equations. Although the term RBFs is used below for Gaussians, some 
of our observations may be relevant also for other types of RBFs. 

1. One of the few methods for constructing efficient function spaces for multi-D quantum problems 
is given by the quasi random distributed Gaussians of [TJ. However, the choice made in [7] for 
the width parameter of an RBF ipi depends only on the potential value at a point, see equation [9l 
while (fi should capture solution properties in a neighborhood of Xj. Such properties correspond 
to the differing nature of the potential function in various regions. Strong repulsions cause 
molecular potentials to rise very steeply when distance coordinates get small. On the other 
hand, the smooth interaction of distant particles causes molecular potentials to be much softer 
for large values of distance coordinates. See for example the Morse potential in figure [H which 
is a standard model for two atom molecular potentials. 

A possible problem in the width choice of [7J can be illustrated by considering two RBFs <pi 
and ipj such that V{xi) = V(xj) are close to the cut off value E max , and such that Xi is in the 
"hard" region of the Morse potential while xj is in the "soft" one. Our assumptions imply that 
<Pi and (pj will have the same large width. In the soft region this is fine, wide basis functions 
are compatible with the smooth character of wavefunctions in this region. However a very wide 
basis function seems unsuitable in the hard region where wavefunction behaviour can change 
from severe oscillation to rapid decay over a short distance. We therefore conclude that RBFs 
centered in hard potential regions should generally be narrower than those centered in soft 
potential regions. This observation is supported by [6] where widths and centers of optimal 
Gaussian RBFs (w.r.t. functional ©) are found for a Morse potential. Comparing centers with 
similiar potential values we see that the width for a center in the hard region is smaller than for 
a center in the soft region (see particularly figure 2 in [6]). Recall that functional ([5|) is different 
from functional ([U]) used for optimizing RBF parameters in [7J . While the reasoning giving ([5]) is 
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clear, we find no similiar justification for ([6]). Particularly, it is not clear why the representation 
of H as a bilinear form, rather than a linear operator, appears in ([6]). 

We suggest that the parameter choice of [7J may be emiliorated by appropriately decreasing the 
widths of RBFs, and increasing their density, in hard potential regions. One way to achieve 
this is to choose nodes and width parameters as in [7J but using a modified potential. Let us 
denote the classically accessible region by C , the "soft" potential region by Q s , and the "hard" 
potential region by Q^. We define V{x) such that V{x) = V(x) for x E £l s n C , and such that 
V(x) has appropriately low values for x G D Sl c - We may choose, for example, the value of 
V{x) in a nearby minimum or in the global one. In the classically inaccessible parts of we set 
V{x) = oo, i.e. we allow no nodes there. The nodes are then chosen similiarly to [7J, but using 
V{x) instead of V{x). Figure Q] illustrates V{x) for a Morse potential. In 1-D and 2-D problems 
this approach should not pose serious difficulty, the "hard" potential regions and nearby minima 
may be identified by inspection of the potential surface. In higher dimensions potential surface 
geometry is not so easy to visualize, yet simple ways may still exist for modifying the choice of 
nodes and width parameters. We defer their discussion to section [5j 

2. In [7], where Gaussian RBFs are used in Galerkin eigenvalue calculations, it is observed that ac- 
curacy improves as the global width parameter c is decreased until unstable calculations appear. 
Therefore cond(S) is used as an indicator for a good choice of c; c is chosen so that cond(S) is in 
the range 10 9 — 10 12 . However, it is desirable that our choice of c will relate to the accuracy of 
our results in a more direct way, and will give a narrower range of possible c values. We therefore 
suggest a different approach for choosing c. Let U be a matrix whose columns give (a few or 
many) approximate eigenfunctions calculated using collocation in the Gaussian RBF basis. Since 
the exact eigenfunctions are orthonormal we suggest to use the "collocation orthogonalization 
error" \\U*SU — I\\ as a performance indicator for the choice of c in the collocation approach. 
To our surprise the values of c minimizing the collocation orthogonalization error were close to 
those minimizing eigenfunction and eigenvalue errors also in the Galerkin approach. Another 
criterion for choosing c based on Hamiltonian trace minimization is suggested in section [5j 

3. The simplicity of RBF collocation is very appealing. But how does collocation compare to the 
traditional Galerkin approach? It is reported in [13], |14j that for their basis functions collocation 
accuracy is similiar to that of Galerkin. They therefore conclude that collocation accuracy is 
similiar to that of Galerkin. However, in all our numerical examples the best Galerkin results 
were better than the best collocation results. 
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Morse oscillator A: potential and eigenfunctions 




- potential 

modified potential 
eigenfunctions 
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Figure 1: The Morse A potential V(x) = 12(1 — e -0 ' 5 ^) 2 is shown here together with the 0th, 10th, and 
23rd eigenfunctions. The graph of each eigenfunction is translated vertically by its own eigenvalue. 
The dotted line is the graph of the corresponding modified potential V(x) defined in equation (|13p . 



Another question is how collocation works together with the quasi randomly distributed Gaus- 
sians of [7], i.e. how does performance compare to Gaussians whose centers are distributed 
deterministically with the same density? To the best of our knowledge this was not previously 
examined in the literature. In our numerical examples collocation performed generally worse 
with the quasi random centers. 



4 Numerical examples and discussion 

The suggestions and observations above were tested on two Hamiltonians which we call Morse A, 
and Morse B. Morse A, from [13], has potential V(x) = 12(1 — e _0 - 5:r ) 2 which is plotted in figure 
[TJ and parameters m = 6, U = 1. This system has 24 bound states. Morse B has the potential 
V(x) = 12(1 — e ~°- 9566:r ) 2 ) anc l a n other parameter values identical to those of Morse A. This system 
has 13 bound states. Analytic formulas for eigenvalues and eigenfunctions were taken from |18j . 

For each Hamiltonian three sets of Gaussian basis functions were tested. All consisted of 48 
Gaussians, the parameter values E max = 12, fj, = 0.02 were used in the construction of all, and all had 
nodes in an interval [L, R]. L is the left classical turning point for E max = 12 with values L ~ —1.3863, 



11 



and L ~ —0.7246, for system A and B respectively. R was chosen arbitrarily (with a few trials), but 
well within the region where the highest bound state has decayed close to zero. The values R = 40, and 
R = 50, were chosen for system A and B respectively (in both cases the right classical turning point 
for E max = 12 is too big for our purposes). The method for choosing nodes and width parameters in 
each set of basis functions follows. 



1. Set 1: Nodes were generated so: x\ = R, Xi+i = X\ — k/\/E max + Jl — V(xi), where k was chosen 
so that X48 = L. This corresponds to the node density of [7J with 7 = 1/2, see equation (J7J). 
Having the nodes, the width parameters were chosen as specified in equation Q. 

2. Set 2: The nodes were generated quasi randomly as specified in [7J with the same density as in 
basis 1. The Sobol sequence software package by W. Putschogl [12] was used. Having the nodes, 
the width parameters were chosen as in basis 1. 

3. Set 3: We regard the left half line x < as the hard potential region for both our Morse 
potentials. Accordingly, set 3 was generated as set 1, but using the following modified potential 
function, which is illustrated in figure [TJ 

V(x) = { V ' . 13 

[ L < x < 

The nodes and width parameters of sets 1A, 2A, 3A, are illustrated in figure [2] (this obvious notation 
distinguishes between basis sets for Hamiltonians A and B). The nodes of set 3A are illustrated in 
figure[3]together with the 23rd eigenfunction, the compatibility between node density and eigenfunction 
oscillation is evident. The node density given by equation © with 7 = 1 was also tested, but results 
were better with 7 = 1/2 as above. 

For each basis set the Hamiltonian was discretized using collocation with the Xi serving as colloca- 
tion nodes, and using Galerkin with analytic calculation of matrix elements. The average eigenvalue 
error 

1 N 



71=0 

and the average L2 norm of eigenfunction errors 



-Y^\E a n -E n \, (14) 



1 x 



n=0 



T Ell^n-^n||L 2 , (15) 



were calculated for different values of the global width parameter c. The superscript a stands for 
"approximate" , and = 23 or 12 for system A or B respectively. Results are shown in figures HI 
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Basis set 1A: nodes/widths 




nodes x. 
widths w 



00 9 a 
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Basis set 
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Basis set 3A: nodes/widths 
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Figure 2: Nodes Xj (dots: •) and width parameters Wi (circles: o) used for the different kinds of 
basis functions. Set 1A: deterministic node choice. Set 2A: quasi random node choice [7j. Set 3A: 
deterministic node choice with density corresponding to the modified potential in equation f)13[) . See 
section [J] for a full explanation of these parameter choices. 



3 set 3A nodes, and 23rd Morse A eigenfunction 




15 20 25 30 35 40 

close up on well region 




Figure 3: Nodes of basis set 3A and 23rd Morse A eigenfunction. Note the compatibility between 
eigenfunction oscillation and node density. The bottom graph shows a close up on the well region. 
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These figures also illustrate the c dependence of the collocation orthogonalization error 

\\U SU — I\\Frobenius , (16) 

with U being the 48 x (iV + 1) matrix of eigenvectors calculated by collocation. Table Q] gives cond(S) 
and cond{$) at the c values minimizing ef for each basis set and discretization method. The only 
exception is in basis 2B were the minimizing c values gave very ill conditioned bases. Instead larger 
c was taken so that cond(S) ~ 10 12 as advocated in [TJ. For the different discretization methods and 
bases (with the c values of table UJ) the logarithms of eigenvalue and eigenfunction errors are plotted 
against their index in figures [5l [7J [8j [H 

Figures [8J and [9] show that basis sets 3 gave orders of magnitude better results compared to basis 
sets 1 and 2, except for a few high eigenfunctions/eigenvalues where accuracy was similiar. This 
indicates that indeed there is room for improvement in the method of [7J for choosing nodes/widths 
in hard potential regions. Moreover, the modification suggested here in basis set 3 can produce highly 
accurate results. 

The results for basis sets 1 and 3 in figures 01 [H (top and bottom rows) indicate that it may 
be possible to choose the global width parameter c based on the collocation orthogonalization error. 
The minimum of e/ for collocation, and for Galerkin, occurs near minima of the collocation "orth 
err". In some cases these minima occur for c values for which cond(S) is not too large, see table 
[TJ This supports our notion that it is desirable to have a criterion for choosing c which is more 
refined than reaching large cond(S). Similiar relations between minima of e/ and of the collocation 
orthogonalization error were not observed in basis sets 2 with the quasi randomly distributed nodes 

m- 

The simplicity of collocation is an important advantage which can save much programing time. 
Moreover, with a careful choice of basis function parameters collocation can give accurate results. 
However, figures [5j [7J illustrate that in all our calculations the Galerkin results were much more 
accurate than those of collocation. 

Examining the bottom row in figures [51 we see that collocation errors were larger with bases 
2A and 2B compared to bases 1A and 2A. That is, we see that larger errors appeared when using 
collocation with quasi random DGBs, compared to DGBs whose centers had the same density but 
were chosen deterministically. 
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5 Conclusions and future directions 



In this work several observations about Gaussian radial basis functions and collocation in Schrodinger 
eigenproblems were made and tested numerically. It may be hoped that the following conclusions are 
relevant not only for choosing good Gaussian bases, but also for other types of localized radial basis 
functions. 

In [7J Garashchuk and Light give one of the few examples of efficient basis functions for multi-D 
Schrodinger equations. An improvement of their approach may be obtained from our observation that 
the Gaussians of [7j are too wide in "hard" potential regions, e.g. where nuclei are close together 
in molecular potentials. This led us to suggest choosing nodes/weights using the method of [7J, but 
with a modified potential function obtained by replacing steeply rising parts by "vertical" walls and 
flattening nearby valley regions. We emphasize that our goal is to produce basis functions which are 
not overly wide and sparse in hard potential regions. Using the modified potential above is just one 
way to achieve this goal, others are suggested in the future directions below. Our numerical tests 
indicate that the changes we propose in the choice of nodes/widths can give a dramatic improvement 
in accuracy. 

In the literature we often encounter the conception that with appropriately wide basis functions 
collocation errors are not larger than those of Galerkin [13], p3], [15]. However, in our numerical 
tests the Galerkin accuracy was better than that of collocation, often orders of magnitude better. The 
difference is probably because we have chosen the optimal global width parameter for each method, 
while in [13], [14] . |15] no special attempt was made to optimize c. Our calculations indicate that 
caution should be exercised before embracing collocation as equally accurate to Galerkin. 

In [TJ, [13], and other places, it is suggested to decrease the global width parameter, making 
basis functions wider, until the basis is nearly badly conditioned. However a more refined method 
for choosing c is desirable. One supporting fact is that the optimal c values in table [T] indeed give 
large cond(S) but not quite as large as expected in the literature. Our numerical tests indicate that 
the collocation orthogonalization error, defined in equation (|16p . may guide us in choosing the global 
width parameter. We found that minima of the collocation orthogonalization error plotted against the 
global width parameter occur near minima of eigenfunction/eigenvalue errors calculated by collocation 
and by the Galerkin approach. 

For the future, it may be profitable to attempt a further refinement of the choice of nodes/widths 
in "hard" potential regions. Suppose we have nodes/widths constructed by some distribution method, 
e.g. that of [7J. Methods should be explored for identifying nodes lying in hard potential regions, 
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for adding more nodes to these regions, and for choosing width parameters with respect to potential 
values in a neighborhood and not just a point. It seems feasible that these goals can be achieved for 
high dimensional systems using only simple local calculations. 

Apart from choosing the centers and relative widths of radial basis functions, a correct scaling of 
their widths through the global width parameter c is critical for good performance. Our suggestion of 
using the collocation orthogonalization error as a criterion for choosing c may be worth pursuing. One 
possibility is to avoid full diagonalization of the collocation Hamiltonian matrix and to measure the 
orthogonalization error only for the few lowest, or highest, eigenvectors. An even simpler method may 
be tried in the Galerkin approach where the "variational principle" holds. In this case the trace of the 
Hamiltonian matrix is minimal for the optimal c (w.r.t. eigenvalue errors). A very simple method for 
choosing a "good" c may possibly be based on this observation. The work required will not be larger 
compared to the method used in [7j based on the condition number of S^ 1 . Some modifications may 
also be examined. For example, if the matrix elements of H are not calculated exactly we may try to 
choose c by minimizing the first few eigenvalues of S -1 !!. 

Recall from section [3] that optimization of functional ([5]) gives nodes/ widths which are consis- 
tent with our observations. So it may be interesting (but probably difficult) to apply the ideas of 
Garashchuk and Light from [7] but with the functional ([5]), rather than the functional ©. That is, to 
attempt identifying a functional form for the densities and width parameters arising from optimization 
of functional © . 

Finally, in addition to further numerical work, a mathematical error analysis for RBFs applied to 
Schrodingers equation should be attempted. 
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Basis set 1A: logfmean eigenvalue error) vs c 



Basis set 1 A: logfmean eigenfu notion error) vs c 




c (global width parameter) 
Basis set 2A: logfmean eigenvalue error) vs c 




c (global width parameter) 
Basis set 3A: log(mean eigenvalue error) vs c 




c (global width parameter) 
Basis set 2A: log(mean eigenfunction error) vs C 



e.f. err galerkin 
e.f. err collocation 
- orth err collocation 



c (global width parameter) 
Basis set 3A: logfmean eigenfunction error) vs c 




c (global width parameter) 



c (global width parameter) 



Figure 4: Error dependence on global width parameter, Morse A: Logarithms of mean eigenvalue 
(left column) and eigenfunction (right column) errors plotted against the global width parameter for 
the different basis sets (rows), and discretization methods, applied to Morse Hamiltonian A. Circles- 
Galerkin, squares-collocation, continuos line- "orth error" which is defined in equation ([T6]) . See table 
[I] A for information on the condition numbers of S and <£. 
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Basis 1 A: log(eigvalue error) vs n 



Basis 2A: log(eigvalue error) vs n 



Basis 3A: log(eigvalue error) vs n 



Basis 1 A: log(eigfunction error) vs n 



o o ° ° ° o o 



Basis 2A: log(eigfunction error) vs n 



Basis 3A: log(eigfunction error) vs n 



Figure 5: Comparing collocation and Galerkin performance, Morse A: Logarithms of eigenvalue (left 
column) and eigenfunction (right column) errors vs. index for each basis set (rows) and each dis- 
cretization method (Circles-Galerkin, squares-collocation) applied to Morse Hamiltonian A. In each 
case results were obtained using the global width parameter value which minimizes the mean eigen- 
function error ej, see table [T]A. Note the low errors in the Galerkin/basis set 3 combination. 
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Basis set 1B: logfmean eigenvalue error) vs c 



Basis set 1 B: logfmean eigenfu notion error) vs c 




c (global width parameter) 
Basis set 2B: logfmean eigenvalue error) vs c 




c (global width parameter) 
Basis set 3B: log(mean eigenvalue error) vs c 



o e.v. err galerkin 
□ e.v. err collocation 
orth err collocation 





c (global width parameter) 
Basis set 2B: log(mean eigenfunction error) vs C 



e.t. err galerkin 
e.f. err collocation 
orth err collocation 




i o o o o 



c (global width parameter) 
Basis set 3B: logfmean eigenfunction error) vs c 



O e.t. err galerkin 
□ e.f. err collocation 
orth err collocation 
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lOOOOOOOOOOOOOOO II 




c (global width parameter) 



c (global width parameter) 



Figure 6: Error dependence on global width parameter, Morse B: Logarithms of mean eigenvalue 
(left column) and eigenfunction (right column) errors plotted against the global width parameter for 
the different basis sets (rows), and discretization methods, applied to Morse Hamiltonian B. Circles- 
Galerkin, squares-collocation, continuos line- "orth error" which is defined in equation (|16p . See table 
Q]B for information on the condition numbers of 5 and <J>. 
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Basis 1B: log(eigvalue error) vs n 



Basis 1B: log(eigfunctior) error) vs n 



Basis 2B: log(eigvalue error) vs n 



Basis 2B: log(eigfunction error) vs n 



Basis 3B: log(eigvalue error) vs n 



Basis 3B: log(eigfunction error) vs n 



Figure 7: Comparing collocation and Galerkin performance, Morse B: Logarithms of eigenvalue (left 
column) and eigenfunction (right column) errors vs. index for each basis set (rows) and each discretiza- 
tion method (Circles-Galerkin, squares-collocation) applied to Morse Hamiltonian B. In each case (with 
the exception of basis 2B) results were obtained using the global width parameter value which mini- 
mizes the mean eigenfunction error e/, see table [I] B. Note the low errors in the Galer kin/basis set 3 
combination. 
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Morse A: Galerkin log(eigvalue error) vs n 



□ basis 2A 
basis 3A 



Morse A: Galerkin log(eigfunction error) vs n 
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□ basis 2A 
basis 3A 



Morse A: Collocation log(eigvalue error) vs n 
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Morse A: Collocation log(eigfunction error) vs n 



O basis 1 A 
□ basis 2A 
basis 3A 



O basis 1 A 
□ basis 2A 
basis 3A 



Figure 8: Comparing the different basis sets, Morse A: Logarithms of eigenvalue (left column) and 
eigenfunction (right column) errors vs. index for each discretization method (rows) and each basis 
set (circles-set 1A, squares-set 2A, diamonds-set 3A) applied to Morse Hamiltonian A. In each case 
results were obtained using the global width parameter value which minimizes the mean eigenfunction 
error ej, see table [T] A. 
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Morse B: Galerkin log{eigvalue error) vs n 



O basis 1B 
□ basis 2B 
basis 3B 



Morse B: Galerkin log(eigfunction error) vs n 



Morse B: Collocation log(eigvalue error) vs n 



o o 



O basis 1 B 
□ basis 2B 
basis 3B 



Morse B: Collocation log(eigfunction error) vs n 



o o 



Figure 9: Comparing the different basis sets, Morse B: Logarithms of eigenvalue (left column) and 
eigenfunction (right column) errors vs. index for each discretization method (rows) and each basis 
set (circles-set IB, squares-set 2B, diamonds-set 3B) applied to Morse Hamiltonian B. In each case 
(with the exception of basis 2B) results were obtained using the global width parameter value which 
minimizes the mean eigenfunction error e/, see table [I] B. 
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basis set 1A 


basis set 2A 


basis set 3A 




Galerkin 


collocation 


Galerkin 


collocation 


Galerkin 


collocation 


c 


0.6 


1.1 


0.8 


1.7 


0.6 


0.9 


cond(Q) 


4.07 x 10 10 


3.99 x 10 3 


6.73 x 10 6 


480 


1.2 x 10 8 


6.53 x 10 4 


cond(S) 


4.4 x 10 9 


8.9 x 10 4 


6.97 x 10 8 


7.3 x 10 4 


4.76 x 10 9 


9.66 x 10 6 




basis set IB 


basis set 2B 


basis set 3B 




Galerkin 


collocation 


Galerkin 


collocation 


Galerkin 


collocation 


c 


1.6 


1.7 


2 


2 


1.4 


1.1 


cond{Q) 


4.12 x 10 7 


6.71 x 10 6 


1.46 x 10 9 


1.46 x 10 9 


5.31 x 10 7 


1.02 x 10 8 


cond(S) 


3.98 x 10 s 


2.21 x 10 s 


1.85 x 10 12 


1.85 x 10 12 


4.28 x 10 9 


1.96 x 10 10 



Table 1: The condition numbers of $ and S are given at the c (global width parameter) values 
minimizing the mean eigenfunction error ej (defined in equation (|15p ) for each basis set and each 
discretization method applied to systems A and B. For basis 2B the non optimal value c = 2 was chosen 
since the associated cond(S) has the maximal value advocated in [7] . The optimal values were c = 0.9 
for Galerkin, with cond(S) = 1.68 x 10 17 , and c = 1.4 for collocation, with cond(S) = 7.79 x 10 14 . 

References 



[1] X.G Wang and T. Carrington Jr. A contracted basis-Lanczos calculation of vibrational levels of 
methane: Solving the Schrddinger equation in nine dimensions, J.Chem.Phys. 119 (2003) 101. 

[2] D.O. Harris, G.G. Engerholm and W.D. Gwinn, Calculation of matrix elements for one- 
dimensional quantum mechanical problems and the application to anharmonics oscillators, 
J.Chem.Phys. 43 (1965) 151; A.S. Dickinson and P.R. Certain, Calculation of matrix elements 
for one- dimensional quantum mechanical problems, J.Chem.Phys. 49 (1968) 4209; R.G. Little- 
john, M. Cargo, T. Carrington Jr., K.A. Mitchell, B. Poirier, A general framework for discrete 
variable representation basis sets, J.Chem.Phys. 116 (2002) 8691. 

[3] J.C. Light and T. Carrington Jr. Discrete- Variable Representations and Their Utilization, Adv. 
in Chem.Phys. 114 (2000) 263, ed. I.Prigogine and S.A. Rice. 

[4] I. Degani and D.J. Tannor, Calculating multidimensional discrete variable representations from 
cubature formulas, J.Phys.Chem. A 110 (2006) 5395; M. Cargo, and R.G. Littlejohn. Tetrahedrally 
invariant discrete variable representation basis on the sphere, J.Chem.Phys. 117 (2002) 59; R. 



23 



Dawes and T. Carrington Jr. A multidimensional discrete variable representation basis obtained 
by simultaneous diagonalization, J.Chem.Phys. 121 (2004) 726. 

[5] LP. Hamilton and J.C. Light, On distributed Gaussian bases for simple model multidimensional 
vibrational problems, J.Chem.Phys. 84 (1986) 306. 

[6] B. Poirier and J.C. Light, Efficient distributed Gaussian basis for rovibrational spectroscopy cal- 
culations, J.Chem.Phys. 113 (2000) 211. 

[7] S. Garashchuk and J.C. Light, Quasirandom distributed Gaussian bases for bound problems, 
J.Chem.Phys. 114 (2001) 3929. 

[8] B.W. Shore, Comparison of matrix methods applied to the radial Schrddinger eigenvalue equa- 
tion: The Morse potential, J.Chem.Phys. 59 (1973) 6450; B.W. Shore, B-spline expansion bases 
for excited states and discretized scattering states, J.Chem.Phys. 63 (1975) 3385; H. Bachau, E. 
Cormier, P. Decleva, J.E. Hansen, and F. Martin, Applications of B- splines in atomic and molec- 
ular physics, Rep.Prog.Phys. 64 (2001) 1815. 

[9] J.K.L. MacDonald, Succesive approximations by the Rayliegh-Ritz variational method, Phs.Rev 43 
(1933) 830. 

[10] G.B. Arfken and H.J. Weber, Mathematical methods for physicists, Harcourt academic press, fifth 
ed. (2001) 817. 

[11] X.G. Hu, T.S. Ho, and H. Rabitz, The collocation method based on a generalized inverse multi- 
quadric basis for bound state problems, Comp.Phys.Comm. 113 (1998) 168. 

[12] W. Putschogl, Interface to Sobol GSL implementation, the matlab central file exchange (2006) 
http://www.mathworks.com/matlabcentral/ 

[13] W. Yang and A.C. Peet, The collocation method for bound solutions of the Schrddinger equation, 
Chem.Phys.Lett. 153 (1988) 98. 

[14] W. Yang and A.C. Peet, The collocation method for calculating vibrational bound states of molec- 
ular systems with application to Ar-HCl, J.Chem.Phys. 90 (1989) 1746. 

[15] W. Yang and A.C. Peet, A method for calculating vibrational bound states: iterative solution of 
the collocation equations constructed from localized basis sets, J.Chem.Phys. 92 (1990) 522. 

[16] L. Collatz, The numerical treatment of differential equations, Springer Verlag, third ed. (1960). 



24 



[17] E.J. Kansa, Multiquadrics - a scattered data approximation scheme with applications to compu- 
tational fluid dynamics - I, Computers Math.Applic. 19 (1990) 127; E.J. Kansa, Multiquadrics - 
a scattered data approximation scheme with applications to computational fluid dynamics - II, 
Computers Math.Applic. 19 (1990) 147. 

[18] M.M. Nieto and L.M. Simmons Jr. Eigenstates, coherent states, and uncertainty products of the 
Morse oscillator, Phys.Rev.A 19 (1979) 438. 

[19] M.D. Buhmann, Radial basis functions, Acta Numerica (2000) 1; R. Schaback and H. Wendland, 
Kernel techniques: From machine learning to meshless methods, Acta Numerica (2006) 543. 



25 



